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Abstract: For the understanding of the origin and propagation of ultra-high energy cosmic rays (UHECR) we developed 
a new approach to simulating UHECRs from an arbitrary number of sources based on Monte Carlo technique. The method 
consists of a combination of three steps. For distant sources we apply commonly accepted parameterizations to calculate 
the contribution to the observed cosmic ray flux. For sources of the local universe we use forward tracking through 
realistic matter distributions and magnetic fields resulting from explicit simulations of large-scale structure formation. 
From the calculations and the forward tracking we generate maps of the probability to observe a particle with a given 
energy from a discrete direction. To account for deflections in the galactic field, these probability maps are transformed by 
matrices calculated from backtracking of antiparticles through field parameterizations. Based on the combined probability 
maps, Monte Carlo production of individual UHECR data is performed which are then used in comparisons with UHECR 
measurements. The simulated UHECR data serves for optimizing the analysis techniques used in UHECR measurements 
as well as for constraining the parameter space of the underlying source and magnetic field models. 
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1 Introduction 

The comparison of observations with model predictions by 
means of Monte Carlo generated data can be a powerful 
tool to improve the knowledge on the origin and propa- 
gation of UHECRs. Currently there are basically two ap- 
proaches to this challenge. 

In the first approach UHECR data are generated by a full 
forward simulation of the UHECR propagation from their 
sources to the observer, accounting for energy losses and 
deflections in magnetic fields at every point of the trajec- 
tory. The low rate of trajectories hitting the observer ren- 
ders the generation of large datasets by this technique a 
computationally intensive task. It becomes more and more 
costly as the source distances increase. Furthermore en- 
ergy losses due to the expansion of the universe and the 
corresponding increase of the photon density are difficult 
to incorporate efficiently. 

A second approach uses parametrizations for the energy 
losses and deflections in magnetic fields and is hence much 
less computationally intensive than the full forward sim- 
ulation. These parametrized simulations however cannot 
easily describe structured extragalactic magnetic fields as 
expected from large-scale structure simulations. Further- 
more the commonly used parametrizations are the result of 
an averaging process over long propagation distances and 
are thus not reliable on short propagation distances. How- 



ever, energy losses due to the expansion of the universe and 
the changes in the photon fields can be easily included. 

In the following we present a combination of both ap- 
proaches for the propagation of protons. For the UHECR 
flux from sources up to a distance of 1 10 Mpc, in the fol- 
lowing called the local universe, we track individual par- 
ticles using matter distributions and magnetic fields from 
constrained simulations of large-scale structure formation. 
The size of this simulation regime is determined by the size 
of the structure simulations at hand. Furthermore a current 
proper distance of 110 Mpc correspond to a red shift of 
z = 0.027 or a scaling factor of i? = 0.97 respectively. 
Neglecting the increased energy losses and only account- 
ing for cosmic time dilation the luminosity of sources in 
the maximum distance is overestimated by at least 3% de- 
pending on the magnetic field strength. 

2 Forward simulation 

The magnetic fields and source distributions used in the for- 
ward simulation are the result from simulations of ACDM 
large scale structure formation |1|. The simulations have 
been constrained using data from the IRAS 1 .2Jy Redshift 
Survey [2| to reproduce the observed distribution of lumi- 
nous matter. These simulations have been combined with 
simulations of the magneto hydrodynamic evolution of a 
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Figure 1: Average deflection strength of the forward simulation of UHECRs from sources in 107 Mpc distance in an 
exemplary local universe and best fit of simulations with PARSEC (left panel). Integrated contribution from sources up 
to distance d to the flux for different energies (right panel). The vertical dashed fine denotes the merging point of the two 
simulation regimes. 



primordial seed field. The resulting magnetic fields are 
shown to be compatible with observations by means of syn- 
thetic Faraday measurements. 

In our simulations UHECRs are injected at point sources 
with a distribution following the distribution of matter in 
the large-scale structure simulation. The trajectories of the 
UHECRs from the source to a spherical observer in the cen- 
ter of the simulation region are calculated using a modified 
version of the CRPropa software |3, 4|. 

The size of the observer dominates the detection rate of 
UHECRs in the forward simulation. Larger observers how- 
ever introduce systematic errors in the deflection angles 
due to the geometric distortion and reduced propagation 
path lengths. To account for theses effects and also en- 
sure a high performance we adapt the observer size to the 
distance of the source such that the systematic error on the 
deflection remains below 2 deg. 

The simulations of the large-scale structure are performed 
in the smoothed particle hydrodynamics formalism. To 
avoid computationally expensive calculations we pre- 
calculate the magnetic field from the smoothed particle 
data on a regular cubic grid with 100 kpc distance between 
nodes. The magnetic field at arbitrary positions is then cal- 
culated using trilinear interpolation between the values at 
the nodes of a cell. The grid spacing of 100 kpc however 
results in a data size of 150 GB for the complete grid. To 
avoid swapping due to the limited memory of typical desk- 
top PCs, the grid is divided into cubes with an edge length 
of 20 Mpc. The resulting memory load of ca. 1 GB for 
the magnetic field can be easily handled by current desktop 
PCs. All UHECRs in one box are first propagated to the 



edges of the cube. Second the next cube is loaded into 
memory and all corresponding UHECRs are processed. 
The processing of the individual cubes is optionally par- 
allelized on several CPUs. To further improve the perfor- 
mance, UHECRs with a propagation distance larger than 
1 .5 times the linear distance between observer and source 
are dropped. Only about 1 per mille of the UHECRs are 
lost owing to this cut [Aj . 

3 Parametrized simulation 

For the propagation in the distant universe we use the PAR- 
SEC software |5|, which uses parametrizations for the en- 
ergy loss and deflections in turbulent fields to calculate a 
map of the probability to observe a particle in a discrete 
direction (pixel) in a given energy bin for extragalactic 
propagation. To account for deflections in galactic mag- 
netic fields, these probability maps are then transformed 
using pre-calculated matrices obtained from backtracking 
of UHECRs through models of the galactic field. Energy 
losses are modelled within continuous energy loss approxi- 
mation implementing the attenuation length from |6| in the 
extragalactic propagation and neglected in our galaxy. 

Assuming source spectra following a power law with spec- 
tral index 7, the probability to observe a particle in pixel j 
with energy Ei in the range < Ei < Eu,i is 

(1) 



32nd International Cosmic Ray Conference, Beijing 201 1 




(b) 



30° 



-120" 



-60° 



60° 



120° 



-60 




Figure 2: Hammer projected probability density maps in galactic coordinates for an exemplary realization with a source 
density of Blue indicates higher, white lower probability on an arbitrary scale, (a) Forward simulation for 

10 EeV. (b) Parametrized simulation for 10 EeV. (c) Weighted combination of map a) and b). (d) Application of a BSS_S 
p9 1 type galactic magnetic field to map c). (e) Probability map from the forward simulation for 50 EeV. (f) Application of 
a BSS_S type galactic magnetic field to map e). 



Here Zk is the redshift of the source at time of emission 
of the particles at source k in current proper distance dk 
with luminosity Lk- Eij and Ei^u denote the energies 
at the source contributing to the observed energies in the 
range Ei^i — £"„ ^ and include energy losses from pho- 
ton interaction and expansion of the universe. The factor 
PEGMF{ctj,k) accounts for the deflection of the UHECR 
in extragalactic fields with aj^k being the angle between 
the direction of pixel j and the direction of source k. 

To describe the deflection in the extragalactic magnetic 
fields we use a Fisher distribution Q to describe the an- 
gular distribution of the UHECRs around the source. The 
Fisher distribution is the normal-distribution on a sphere 
and has the functional form 



We use K — with 
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with concentration parameter k indicating the width of the 
distribution and a the angular distance from the center of 
the distribution. For small concentration parameters k the 
Fisher distribution converges to the isotropic distribution. 
For large k the Fisher distribution converges to a Rayleigh 
distribution with width a and k — l/a^. 
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taken from ||8l to describe the deflection power in mag- 
netic fields. Here Dk is the distance of source fc, A the 
coherence length, and B the strength of the magnetic field. 
E denotes the energy of the cosmic ray. The elongation 
of the propagation path due to deflection in the magnetic 
fields is parametrized based on the work of Achterberg et 
al. im, but has been modified to include energy losses in 
the parametrization IS) . 

4 Combination of simulation regimes 

For a consistent model, the deflection strength of the 
parametrized simulation has to be scaled to match the mean 
deflection in the forward propagation. The left panel of fig- 
ure[T]shows the mean deflection resulting from the forward 
simulation of 5000 UHECRs from isotropically distributed 
starting points in 107 Mpc distance in the local universe 
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Figure 3: Energy spectrum resulting from the exem- 
plary simulation in comparison with observations from the 
HIRES experiment lilOJ and the Pierre- Auger Observatory 

m. 

and the best fit of the mean deflection strength with PAR- 
SEC. The best fit is achieved for B^/X = 0.94 nG^/W^ 
in equation |3] 

Both simulation regimes are combined by means of sum- 
ming up maps of the probability to observe a particle in 
a given energy range from a discrete direction. While the 
probability maps of the contribution from the distant uni- 
verse are calculated explicitly with PARSEC, the probabil- 
ity maps from the local universe are created from the for- 
ward tracking data. The maps are weighted so that the total 
contribution from the local universe matches the expected 
contribution from continuous sources in the same volume 
for the lowest energy bin. This weight is calculated using 
PARSEC with the matched strength for the magnetic field. 
The integrated relative contributions over the distance are 
shown in the right panel of figure[T]for three different ener- 
gies . 

To account for deflections in the galactic field, the result- 
ing probability maps are transformed using the matrix tech- 
nique provided by PARSEC. From the final probability 
maps simulated data sets with large numbers of UHECRs 
are then generated with low computational effort. 



verse and more than 4.4 million sources up to the maximum 
simulation distance of 1500 Mpc. The source spectra have 
been set to 7 — —2.7. For an energy of 10 EeV the ratio 
of the flux from the distant universe to the local universe 
Nnear / Ntot — 0-23. For energies above 100 EeV all parti- 
cles originate from the local universe. 

In figure [3] we show the resulting energy spectrum from 
the exemplary realization calculated using 100 energy bins 
between 10^^'^ eV and 10^° eV. For comparison the ob- 
served energy spectra of the HiRes experiment [ 10 1 and the 
Pierre- Auger Observatory IfTTI are also shown. The simu- 
lated spectrum is normalized to match the observation of 
the Pierre Auger Observatory at 10 EeV. 

With the technique presented in this contribution we have 
developed a promising mechanism for the extensive pro- 
duction of Monte Carlo simulations of UHECRs suitable 
for studies of anisotropy and cosmic magnetic fields. The 
simulation includes effects of structured extragalactic and 
galactic magnetic fields, as well as energy losses and arbi- 
trary source configurations. 
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5 Results and Conclusion 

In figure |2] we show exemplary probability maps from 
the forward propagation in the local universe and the 
parametrized simulation in the distant universe together 
with the combined maps at two energies. The source den- 
sity in this exemplary simulation is p = lO^'^-^ Mpc^'^, 
resulting in a total number of 1763 sources in the local uni- 



